################################################################################################################################################################################################################################################
# The Political Consequences of Ethnically Targeted Incarceration: Evidence from Japanese-American Internment During WWII
# Replication Code: SI Simulated Locations from WRA Data
################################################################################################################################################################################################################################################

################################################################################################################################################################################################################################################
# Setup
################################################################################################################################################################################################################################################
# Packages
  require(haven)
  require(data.table)
  require(stargazer)
  require(ggplot2)
  require(nnet)
  require(ggjoy)
  require(RColorBrewer)

# Functions
  splitit <- function(x,splitchar,n) sapply(strsplit(as.character(x), splitchar), "[", n)
  getmode <- function(v) {
    uniqv <- unique(v)
    uniqv[which.max(tabulate(match(v, uniqv)))]
  }
  as.numeric.factor <- function(x) {as.numeric(levels(x))[x]}
  
  cl  <- function(dat, fm, cluster){
    require(sandwich, quietly = TRUE)
    require(lmtest, quietly = TRUE)
    M <- length(unique(cluster))
    N <- length(cluster)
    K <- fm$rank
    dfc <- (M / (M - 1)) * ((N - 1) / (N - K))
    uj  <- apply(estfun(fm), 2, function(x) tapply(x, cluster, sum));
    vcovCL <- dfc * sandwich(fm, meat = crossprod(uj) / N)
    coeftest(fm, vcovCL)}
  
  set.seed(14604)

################################################################################################################################################################################################################################################
# Directories
################################################################################################################################################################################################################################################
# This sets your working directory to the location of the source file. Note that this command is specific to RStudio. You'll need to change
# this if you're working directly in R or some other IDE (see https://stackoverflow.com/questions/13672720/r-command-for-setting-working-directory-to-source-file-location-in-rstudio)
  setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
  
################################################################################################################################################################################################################################################
# Load Data
################################################################################################################################################################################################################################################
# JARP Survey
  full <- readRDS("jarpfull_plus.RDS")

# WRA Data on Interned Population  
  wra <- read_dta("internment_raw.dta")
  wra <- data.frame(wra)

################################################################################################################################################################################################################################################
# Merge in Camp Information for Treatments
################################################################################################################################################################################################################################################
  camps <- sort(unique(full$camp[full$camp != "No Data" & full$camp != "Not evacuated or interned" & is.na(full$camp) == 0]))
  
  # Just looked counties up manually for these in order. Some places are in multiple counties, so the county listed is the one that contains most of the camp
  # Key to major camps: 
  # Gila River = Rivers, AZ in Pinal County, AZ
  # Granada = Amache, CO in Prowers County, CO
  # Heart Mountain = Heart Mountain, Wyoming in Park County
  # Jerome = Denson, AR mostly in Drew but partly in Chicot county
  # Manzanar = Manzanar, CA in Inyo County
  # Minidoka = Hunt, ID in Jerome County, Idaho
  # Poston = Poston, AZ in Yuma County
  # Rohwer = McGehee, AR in Desha County Arkansas
  # Topaz = Topaz, UT in Utah
  # Newell, CA is Tule Lake, which is mostly Modoc but partly in Siskiyou county
  # Camp Livingston, LA is on the Rapides Parish and Grant Parish Lines
  
  # Merge in County Names for Various Camps
  camps <- data.frame(camps = camps, V1 = cbind(splitit(camps, ",", 1)))
  camps$counties <- c("Prowers", "Burleigh", "Zavala", "Drew", "Comanche", "Park", "Jerome", "Rapides", "Hidalgo", "Inyo", "Desha", "Missoula", "Modoc", "Yuma", "Pinal", "Santa Fe", "Millard")
  camps$stateabb <- splitit(camps$camps, ", ", 2)
  camps <- merge(camps, data.frame(cbind(state.abb, state.name)), by.x = "stateabb", by.y = "state.abb")
  camps <- camps[, c("camps", "counties", "state.name")]
  colnames(camps) <- c("camp", "county", "state")
  camps$county <- toupper(camps$county) ; camps$state <- toupper(camps$state)
  
  # Identify Additional Camp Treatments
  # Watch Towers, Military Police Buildings, Strikes (see Confinement and Ethnicity book in References)
  # towers = watch towers in camp; bldgs = military buildings on premises; strikes = strikes or work stoppages;
  # force = camp personnel use force or violence against internees; violence = violence among internees or between internees and civilian population;
  # peakpop = peak population (Table 3.2 in Confinement and Ethnicity)
  # Other sources: http://newmexicohistory.org/places/lordsburg-internment-pow-camp (Lordsburg), http://encyclopedia.densho.org/Fort_Sill_(detention_facility)/
  # 
  severity <- data.frame(camp = c("Rivers, AZ", "Amache, CO", "Heart Mountain, WY", "Denson, AR", "Manzanar, CA", "Hunt, ID",
                                  "Poston, AZ", "McGehee, AR", "Topaz, UT", "Newell, CA", "Fort Sill, OK", "Lordsburg, NM", "Bismarck, ND",
                                  "Crystal City, TX", "Livingston, LA"),
                         campname = c("Gila River", "Amache", "Heart Mountain", "Jerome", "Manzanar", "Minidoka",
                                      "Poston", "Rohwer", "Topaz", "Tule Lake", "Fort Sill", "Lordsburg", "Fort Lincoln",
                                      "Crystal City", "Livingston"),
                         towers = c(1, 6, 9, 7, 8, 8, 0, 8, 7, 19, NA, NA, NA, NA, NA),
                         bldgs = c(15, 15, 19, 12, 13, 15, 12, 12, NA, 63, NA, NA, NA, NA, NA),
                         strikes = c(0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1),
                         force = c(0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0),
                         violence = c(0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),
                         peakpop = c(13348, 7318, 10767, 8497, 10046, 9397, 17814, 8475, 8130, 18789, 707, 2500, 1518, 4000, 1123))
  
  severity$militarism <- (severity$towers) / (severity$peakpop / 1000)
  
  severity$stateview <- ifelse(severity$strikes == 1 | severity$force == 1, 1, 0)
  
  severity$any <- ifelse(severity$strikes == 1 | severity$force == 1 | severity$strikes == 1, 1, 0)
  
  camps <- merge(camps, severity, by = "camp")
  
  full <- merge(full, camps, by = "camp", all.x = T)
  
  
  # Define Treatment Status    
  full$anyfam <- ifelse(full$internedinfam > 0, 1, 0)
  
  full$treatmentstatus <- ifelse(full$anyfam == 1 & full$interned.numeric == "Yes", "self.fam",
                                 ifelse(full$anyfam == 0 & full$interned.numeric == "Yes", "self.only",
                                        ifelse(full$anyfam == 1 & full$interned.numeric == "No", "fam.only",
                                               ifelse(full$anyfam == 0 & full$interned.numeric == "No", "control", NA))))
  
  
  full$exposure <- ifelse(full$treatmentstatus %in% c("self.only", "self.fam"), "self.interned", 
                          ifelse(full$treatmentstatus %in% c("fam.only"), "fam.only", "control"))
  
################################################################################################################################################################################################################################################
# Assign Hypothetical Locations Based on Empirical Locations in WRA
################################################################################################################################################################################################################################################
# Get Empirical Probabilities of Particular County Codes by State
  # WRA State Codes
    wra.state <- read.csv("wra_state.csv", stringsAsFactors = F)
  
  # WRA Address Codes
    wra.addy <- read.csv("wra_state_county.csv", stringsAsFactors = F)  
    colnames(wra.addy) <- c("address", "location", "comments")
    wra.addy <- wra.addy[, !(colnames(wra.addy) %in% "comments"),]
    
    wra <- merge(wra, wra.addy, by = "address", all.x = T)
  
  # Identify Pre-Internment States in WRA Data
    wra$origin_state <- substr(wra$address, 1, 2)
    wra <- merge(wra, wra.state, by.x = "origin_state", by.y = "state", all.x = T)
    wra$state_name <- splitit(wra$state.desc, " - ", 2)
    wra$state_name[grepl("County", wra$state_name) == 1 | grepl("Hawaii", wra$state_name) == 1] <- "Hawaii"
    
  # Get Counts of Internees by County Within State
    wra.locs <- data.table(wra[!is.na(wra$address) & !is.na(wra$state_name),])
    wra.locs <- wra.locs[ , list(counts = .N), by = list(state_name, address)]
    wra.locs <- wra.locs[ , list(props = counts / sum(counts), address = address), by = list(state_name)]
    
  # Assign correlated numeric codes instead of addresses; this will let us regress numerically without
  # too much loss of power - I'm just recycling the JARP's codes here to match our previous analysis closely
    wrastates <- data.frame(state_name = unique(wra.locs$state_name), statecode = NA)
    wrastates$statecode <- c("91", "92", "93", "86", "82", "83", "84", "85", "81", "87", "61", "71", "74", "35",
                             "31", "33", "44", "45", "46", "58", "53", "22", "23", "99", "98", "01")
    
    wra.locs <- merge(wra.locs, wrastates, by = "state_name")
    wra.locs$addcode <- rep(NA, nrow(wra.locs))
    
    for(i in 1:nrow(wrastates)){
      wra.locs$addcode[wra.locs$state_name == wrastates$state_name[i]] <- 1:length(wra.locs$addcode[wra.locs$state_name == wrastates$state_name[i]])
    }
    
    wra.locs$addcode <- as.character(wra.locs$addcode)
    wra.locs$address2 <- ifelse(nchar(wra.locs$addcode) == 1, paste(wra.locs$statecode, "00", wra.locs$addcode, sep = ""),
                                ifelse(nchar(wra.locs$addcode) == 2, paste(wra.locs$statecode, "0", wra.locs$addcode, sep = ""),
                                       ifelse(nchar(wra.locs$addcode) == 3, paste(wra.locs$statecode, wra.locs$addcode, sep = ""), NA)))
    
################################################################################################################################################################################################################################################
# Recode Outcomes & Covariates
################################################################################################################################################################################################################################################
# Demographics
  # Age
    full$age <- as.numeric(full$age)

  # Generate indicator for pre-internment location across generations
    full$i.precamploc <- as.character(full$i.precamploc) ; full$n.precamp <- as.character(full$n.precamp) ; full$s.precamp <- as.character(full$s.precamp)
    
    full$precamp <- ifelse(full$generation == "Issei", full$i.precamploc,
                           ifelse(full$generation == "Nisei", full$n.precamp,
                                  ifelse(full$generation == "Sansei", full$s.precamp, NA)))
  
  # Generate indicator for farm vs. nonfarm occupation. For Issei, this is going to be the head of household's occupation from 1932-1941;
  # for Nisei it will be the respondent's own primary occupation from 1932-1941
  # Convert to numeric  
    full$i.pthhjob <- as.numeric.factor(full$i.pthhjob) ; full$n.ptjob <- as.numeric.factor(full$n.ptjob)
  
  # Generate single indicator because these are both coded via note 3
    full$ptjob <- ifelse(full$generation == "Issei", full$i.pthhjob,
                         ifelse(full$generation == "Nisei", full$n.ptjob, NA))
    
  # Flag farm occupations
    full$ptfarm <- ifelse(is.na(full$ptjob) == 1, NA,
                          ifelse(full$ptjob %in% c(222, 901, 905, 981, 982, 983, 984, 985, 986, 987, 988, 989), 1, 0))
    
# Identify Precamp State in JARP
    locations <- read.csv("postcamp_location.csv", stringsAsFactors = F)
    locations$precamp.state.code <- substr(as.character(locations$postcamp.code), 1, 2)
    locations <- locations[ , c("postcamp.state", "precamp.state.code")]
    locations <- unique(locations)
    colnames(locations)[1] <- "precamp.state.name"
    
    full$precamp.state <- substr(full$precamp, 1, 2)
    
    full <- merge(full, locations, by.x = "precamp.state", by.y = "precamp.state.code", all.x = T)  
    
################################################################################################################################################################################################################################################
# Figure Out How Many Families Actually Live Together
################################################################################################################################################################################################################################################
# Splitting out direct and indirect effects
  fams <- data.table(full[full$exposure %in% c("self.interned", "fam.only"),])
  fams <- fams[ , list(famlocs = length(unique(precamp))), by = list(respid)]
  setkey(fams, famlocs)
  
################################################################################################################################################################################################################################################
# Assign Draws
################################################################################################################################################################################################################################################
# Recode Variables for Regressions
  # Interest
    interest <- unique(c(full$i.interestus, full$n.interestus, full$s.interestus))
    interestcodes <- rep(NA, length(interest))
    
    interest <- data.frame(cbind(interest, interestcodes))
    interest$interestcodes <- ifelse(interest$interest == "No" | interest$interest == "No interest at all", 0,
                                     ifelse(interest$interest == "Only a little" | interest$interest == "Yes", 1,
                                            ifelse(interest$interest == "A fair amount", 2,
                                                   ifelse(interest$interest == "A great deal", 3, NA))))
    
    interest$interest <- as.character(interest$interest)
    
    full$ptxinterest <- ifelse(full$generation == "Issei", full$i.interestus,
                               ifelse(full$generation == "Nisei", full$n.interestus, 
                                      ifelse(full$generation == "Sansei", full$s.interestus, NA)))
    
    full <- merge(full, interest, by.x = "ptxinterest", by.y = "interest")
    
  # Faith in Government    
  # This question is phrased "Most people in government are not really interested in the problems of the average man". So agreeing here means you don't think much of government.
    full$govhelp <- ifelse(full$generation == "Nisei", full$n.govhelp,
                           ifelse(full$generation == "Sansei", full$s.govhelp, NA))
    
    # Numeric coding
    full$govhelp <- ifelse(full$govhelp == "Agree", 1,
                           ifelse(full$govhelp == "Disagree", 0, NA))
    
  # Political Advice  
    # Only asked of Sansei and Nisei
    full$poladvc <- ifelse(full$generation == "Nisei", full$n.poladv,
                           ifelse(full$generation == "Sansei", full$s.poladv, NA))
    
    full$poladvc <- ifelse(full$poladvc == "Yes", 1,
                           ifelse(full$poladvc == "No", 0, NA))  
    
  # Leadership  
    # Generate Combined Leadership Questions
    full$leader <- ifelse(full$generation == "Nisei", full$n.leader, ifelse(full$generation == "Sansei", full$s.leader, NA))
    full$leader <- ifelse(full$leader == "Protest", -1, ifelse(full$leader == "Orderly and comfortable", 1, ifelse(full$leader == "Can't generalize, both, neither" | full$leader == "Can't generalize; both; neither", 0, NA)))
    
# Set # Simulations
  sims <- 10000
  
# Holder vectors for coefficients  
  interest.violence.coefs <- rep(NA, sims)
  interest.violence.inter.coefs <- rep(NA, sims)
  interest.force.coefs <- rep(NA, sims)
  interest.force.inter.coefs <- rep(NA, sims)
  interest.strikes.coefs <- rep(NA, sims)
  interest.strikes.inter.coefs <- rep(NA, sims)
  interest.militarism.coefs <- rep(NA, sims)
  interest.militarism.inter.coefs <- rep(NA, sims)
  advice.violence.coefs <- rep(NA, sims)
  advice.violence.inter.coefs <- rep(NA, sims)
  advice.force.coefs <- rep(NA, sims)
  advice.force.inter.coefs <- rep(NA, sims)
  advice.strikes.coefs <- rep(NA, sims)
  advice.strikes.inter.coefs <- rep(NA, sims)
  advice.militarism.coefs <- rep(NA, sims)
  advice.militarism.inter.coefs <- rep(NA, sims)
  faith.violence.coefs <- rep(NA, sims)
  faith.violence.inter.coefs <- rep(NA, sims)
  faith.force.coefs <- rep(NA, sims)
  faith.force.inter.coefs <- rep(NA, sims)
  faith.strikes.coefs <- rep(NA, sims)
  faith.strikes.inter.coefs <- rep(NA, sims)
  faith.militarism.coefs <- rep(NA, sims)
  faith.militarism.inter.coefs <- rep(NA, sims)
  leadership.violence.coefs <- rep(NA, sims)
  leadership.violence.inter.coefs <- rep(NA, sims)
  leadership.force.coefs <- rep(NA, sims)
  leadership.force.inter.coefs <- rep(NA, sims)
  leadership.strikes.coefs <- rep(NA, sims)
  leadership.strikes.inter.coefs <- rep(NA, sims)
  leadership.militarism.coefs <- rep(NA, sims)
  leadership.militarism.inter.coefs <- rep(NA, sims)
  
# Identify unique precamp locations list
  precamp.locs <- unique(full$precamp.state.name[full$exposure %in% c("self.interned", "fam.only")])
  precamp.locs <- precamp.locs[!is.na(precamp.locs)]

  
# Initialize full simulation loop  
  for(i in 1:sims){
    # Empty column of simulated locations
    full$wra.precamploc <- rep(NA, nrow(full))
    
    # Loop through locations
      for(j in 1:length(precamp.locs)){
        # Grab everyone in the data in a particular precamp.stat
        temp <- full[which(full$precamp.state.name == precamp.locs[j]),]
        
        # Loop through families to assign locations
        fams <- unique(temp$respid)
        for(z in 1:length(fams)){
          # Draw a location based on the family's state
          if(length(wra.locs$props[wra.locs$state_name == precamp.locs[j]]) == 0){
            temp$wra.precamploc[temp$respid == fams[z]] <- NA
          } else if(length(wra.locs$props[wra.locs$state_name == precamp.locs[j]]) == 1){
            temp$wra.precamploc[temp$respid == fams[z]] <- wra.locs$address2[wra.locs$state_name == precamp.locs[j]]
          } else {
          draw <- rmultinom(1, 1, wra.locs$props[wra.locs$state_name == precamp.locs[j]])
          address <- wra.locs$address2[wra.locs$state_name == precamp.locs[j]][which(draw == 1)]
          temp$wra.precamploc[temp$respid == fams[z]] <- address
          }
        }
        
        
      # Put specific address information back into full data
        full$wra.precamploc[which(full$precamp.state.name == precamp.locs[j])] <- temp$wra.precamploc
        
      }
    
    # Run Models on Modified Data and Store Coefficients
      # Interest: Violence
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("interestcodes", "violence", "age", "gender", "wra.precamploc", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        interest.violence <- lm(interestcodes ~ violence * exposure + wra.precamploc + gender + age + generation, data = comp)  
        interest.violence.coefs[i] <- summary(interest.violence)$coefficients["violence","Estimate"]
        interest.violence.inter.coefs[i] <- summary(interest.violence)$coefficients["violence:exposureself.interned","Estimate"]
        
      # Interest: Force  
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("interestcodes", "force", "age", "gender", "wra.precamploc", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        interest.force <- lm(interestcodes ~ force * exposure + wra.precamploc + gender + age + generation, data = comp)  
        interest.force.coefs[i] <- summary(interest.force)$coefficients["force","Estimate"]
        interest.force.inter.coefs[i] <- summary(interest.force)$coefficients["force:exposureself.interned","Estimate"]
        
      # Interest: Strikes
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("interestcodes", "strikes", "age", "gender", "wra.precamploc", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        interest.strikes <- lm(interestcodes ~ strikes * exposure + wra.precamploc + gender + age + generation, data = comp)  
        interest.strikes.coefs[i] <- summary(interest.strikes)$coefficients["strikes","Estimate"]
        interest.strikes.inter.coefs[i] <- summary(interest.strikes)$coefficients["strikes:exposureself.interned","Estimate"]
        
      # Interest: Militarism
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("interestcodes", "militarism", "age", "gender", "wra.precamploc", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        interest.militarism <- lm(interestcodes ~ militarism * exposure + wra.precamploc + gender + age + generation, data = comp)  
        interest.militarism.coefs[i] <- summary(interest.militarism)$coefficients["militarism","Estimate"]
        interest.militarism.inter.coefs[i] <- summary(interest.militarism)$coefficients["militarism:exposureself.interned","Estimate"]
        
      # Faith: Violence  
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("govhelp", "violence", "wra.precamploc", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        faith.violence <- lm(govhelp ~ violence * exposure + wra.precamploc + age + gender + generation,  data = comp)
        faith.violence.coefs[i] <- summary(faith.violence)$coefficients["violence","Estimate"]
        faith.violence.inter.coefs[i] <- summary(faith.violence)$coefficients["violence:exposureself.interned","Estimate"]
      
      # Faith: Force  
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("govhelp", "force", "wra.precamploc", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        faith.force <- lm(govhelp ~ force * exposure + wra.precamploc + age + gender + generation,  data = comp)
        faith.force.coefs[i] <- summary(faith.force)$coefficients["force","Estimate"]
        faith.force.inter.coefs[i] <- summary(faith.force)$coefficients["force:exposureself.interned","Estimate"]
      
      # Faith: Strikes
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("govhelp", "strikes", "wra.precamploc", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        faith.strikes <- lm(govhelp ~ strikes * exposure + wra.precamploc + age + gender + generation,  data = comp)
        faith.strikes.coefs[i] <- summary(faith.strikes)$coefficients["strikes","Estimate"]
        faith.strikes.inter.coefs[i] <- summary(faith.strikes)$coefficients["strikes:exposureself.interned","Estimate"]
        
      # Faith: Militarism
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("govhelp", "militarism", "wra.precamploc", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        faith.militarism <- lm(govhelp ~ militarism * exposure + wra.precamploc + age + gender + generation,  data = comp)
        faith.militarism.coefs[i] <- summary(faith.militarism)$coefficients["militarism","Estimate"]
        faith.militarism.inter.coefs[i] <- summary(faith.militarism)$coefficients["militarism:exposureself.interned","Estimate"]
        
      # Political Advice: Violence  
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("poladvc", "violence", "wra.precamploc", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        polad.violence <- lm(poladvc ~ violence * exposure + wra.precamploc + age + gender + generation,  data = comp)
        advice.violence.coefs[i] <- summary(polad.violence)$coefficients["violence", "Estimate"]
        advice.violence.inter.coefs[i] <- summary(polad.violence)$coefficients["violence:exposureself.interned", "Estimate"]
        
      # Political Advice: Force  
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("poladvc", "force", "wra.precamploc", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        polad.force <- lm(poladvc ~ force * exposure + wra.precamploc + age + gender + generation,  data = comp)
        advice.force.coefs[i] <- summary(polad.force)$coefficients["force", "Estimate"]  
        advice.force.inter.coefs[i] <- summary(polad.force)$coefficients["force:exposureself.interned", "Estimate"]  
        
      # Political Advice: Strikes
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("poladvc", "strikes", "wra.precamploc", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        polad.strikes <- lm(poladvc ~ strikes * exposure + wra.precamploc + age + gender + generation,  data = comp)
        advice.strikes.coefs[i] <- summary(polad.strikes)$coefficients["strikes", "Estimate"]  
        advice.strikes.inter.coefs[i] <- summary(polad.strikes)$coefficients["strikes:exposureself.interned", "Estimate"]  
        
      # Political Advice: Militarism
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("poladvc", "militarism", "wra.precamploc", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        polad.militarism <- lm(poladvc ~ militarism * exposure + wra.precamploc + age + gender + generation,  data = comp)
        advice.militarism.coefs[i] <- summary(polad.militarism)$coefficients["militarism", "Estimate"]  
        advice.militarism.inter.coefs[i] <- summary(polad.militarism)$coefficients["militarism:exposureself.interned", "Estimate"]  
        
      # Leadership: Violence  
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("leader", "violence", "wra.precamploc", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        leader.violence <- lm(leader ~ violence * exposure + wra.precamploc + gender + age + generation, data = comp)  
        leadership.violence.coefs[i] <- summary(leader.violence)$coefficients["violence", "Estimate"]
        leadership.violence.inter.coefs[i] <- summary(leader.violence)$coefficients["violence:exposureself.interned", "Estimate"]
        
      # Leadership: Force  
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("leader", "force", "wra.precamploc", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        leader.force <- lm(leader ~ force * exposure + wra.precamploc + gender + age + generation, data = comp)  
        leadership.force.coefs[i] <- summary(leader.force)$coefficients["force", "Estimate"]
        leadership.force.inter.coefs[i] <- summary(leader.force)$coefficients["force:exposureself.interned", "Estimate"]
        
      # Leadership: Strikes 
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("leader", "strikes", "wra.precamploc", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        leader.strikes <- lm(leader ~ strikes * exposure + wra.precamploc + gender + age + generation, data = comp)  
        leadership.strikes.coefs[i] <- summary(leader.strikes)$coefficients["strikes", "Estimate"]
        leadership.strikes.inter.coefs[i] <- summary(leader.strikes)$coefficients["strikes:exposureself.interned", "Estimate"]
  
      # Leadership: Militarism
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("leader", "militarism", "wra.precamploc", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        leader.militarism <- lm(leader ~ militarism * exposure + wra.precamploc + gender + age + generation, data = comp)  
        leadership.militarism.coefs[i] <- summary(leader.militarism)$coefficients["militarism", "Estimate"]
        leadership.militarism.inter.coefs[i] <- summary(leader.militarism)$coefficients["militarism:exposureself.interned", "Estimate"]
        
        if(i %% 100 == 0){print(i)}
      }
  
  
      wra.robust <- data.frame(results = c(interest.violence.coefs, interest.force.coefs, interest.strikes.coefs, interest.militarism.coefs, interest.violence.coefs + interest.violence.inter.coefs, interest.force.coefs + interest.force.inter.coefs, interest.strikes.coefs + interest.strikes.inter.coefs, interest.militarism.coefs + interest.militarism.inter.coefs,
                                           faith.violence.coefs, faith.force.coefs, faith.strikes.coefs, faith.militarism.coefs, faith.violence.coefs + faith.violence.inter.coefs, faith.force.coefs + faith.force.inter.coefs, faith.strikes.coefs + faith.strikes.inter.coefs, faith.militarism.coefs + faith.militarism.inter.coefs,
                                           leadership.violence.coefs, leadership.force.coefs, leadership.strikes.coefs, leadership.militarism.coefs, leadership.violence.coefs + leadership.violence.inter.coefs, leadership.force.coefs + leadership.force.inter.coefs, leadership.strikes.coefs + leadership.strikes.inter.coefs, leadership.militarism.coefs + leadership.militarism.inter.coefs,
                                           advice.violence.coefs, advice.force.coefs, advice.strikes.coefs, advice.militarism.coefs, advice.violence.coefs + advice.violence.inter.coefs, advice.force.coefs + advice.force.inter.coefs, advice.strikes.coefs + advice.strikes.inter.coefs, advice.militarism.coefs + advice.militarism.inter.coefs),
                               labels = rep(c(rep("Violence", sims), rep("Force", sims), rep("Strikes", sims), rep("Militarism", sims)), 8),
                               groups = c(rep(c(rep("Family-Only Exposure", 40000), rep("Direct Exposure", 40000)), 4)),
                               outcome = c(rep("Political Interest", 80000), rep("Political Distrust", 80000), rep("Leadership Approach", 80000), rep("Political Advice", 80000)))
      
      wra.robust.plot <- data.table(wra.robust)
      
      wra.robust.plot <- data.frame(wra.robust.plot[, list(result = mean(results), lb95 = quantile(results, 0.025), ub95 = quantile(results, 0.975),
                             lb84 = quantile(results, 0.08), ub84 = quantile(results, 0.092)), by = list(labels, groups, outcome)])
      
      wra.robust.plot$groups <- factor(wra.robust.plot$groups, levels = c( "Family-Only Exposure", "Direct Exposure"), ordered = T)
      
      pdf("si_fig12a.pdf")
      ggplot(data = wra.robust.plot[wra.robust.plot$labels == "Strikes",], aes(x = outcome, y = result, color = groups)) + 
        geom_hline(yintercept = 0, alpha = .5) + 
        geom_pointrange(data = wra.robust.plot[wra.robust.plot$labels == "Strikes",], aes(ymin = lb95, ymax = ub95), size = .5, position = position_dodge(.5)) + 
        geom_pointrange(data = wra.robust.plot[wra.robust.plot$labels == "Strikes",], aes(ymin = lb84, ymax = ub84), size = .5, position = position_dodge(.5)) + 
        theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) + 
        labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Witnessing Demonstrations While Interned", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
      dev.off()
      
      pdf("si_fig12b.pdf")
      ggplot(data = wra.robust.plot[wra.robust.plot$labels == "Violence",], aes(x = outcome, y = result, color = groups)) + 
        geom_hline(yintercept = 0, alpha = .5) + 
        geom_pointrange(data = wra.robust.plot[wra.robust.plot$labels == "Violence",], aes(ymin = lb95, ymax = ub95), size = .5, position = position_dodge(.5)) + 
        geom_pointrange(data = wra.robust.plot[wra.robust.plot$labels == "Violence",], aes(ymin = lb84, ymax = ub84), size = .5, position = position_dodge(.5)) + 
        theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) + 
        labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Witnessing Violence While Interned", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
      dev.off()
      
      pdf("si_fig13a.pdf")
      ggplot(data = wra.robust.plot[wra.robust.plot$labels == "Force",], aes(x = outcome, y = result, color = groups)) + 
        geom_hline(yintercept = 0, alpha = .5) + 
        geom_pointrange(data = wra.robust.plot[wra.robust.plot$labels == "Force",], aes(ymin = lb95, ymax = ub95), size = .5, position = position_dodge(.5)) + 
        geom_pointrange(data = wra.robust.plot[wra.robust.plot$labels == "Force",], aes(ymin = lb84, ymax = ub84), size = .5, position = position_dodge(.5)) + 
        theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) + 
        labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Witnessing Force While Interned", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
      dev.off()
      
      pdf("si_fig13b.pdf")
      ggplot(data = wra.robust.plot[wra.robust.plot$labels == "Militarism",], aes(x = outcome, y = result, color = groups)) + 
        geom_hline(yintercept = 0, alpha = .5) + 
        geom_pointrange(data = wra.robust.plot[wra.robust.plot$labels == "Militarism",], aes(ymin = lb95, ymax = ub95), size = .5, position = position_dodge(.5)) + 
        geom_pointrange(data = wra.robust.plot[wra.robust.plot$labels == "Militarism",], aes(ymin = lb84, ymax = ub84), size = .5, position = position_dodge(.5)) + 
        theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) + 
        labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Militarized Camp Environment", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
      dev.off()
      